Restoring the density-gradient expansion for exchange in solids and surfaces 
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Successful modern generalized gradient approximations (GGA's) are biased toward atomic ener- 
gies. Restoration of the first-principles gradient expansion for exchange over a wide range of density 
gradients eliminates this bias. We introduce PBEsol, a revised Perdew-Burke-Ernzerhof GGA that 
improves equilibrium properties of densely-packed solids and their surfaces. 
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Ground-state Kohn-Sham density functional theory 
(DPT) [T] has been hugely successful for electronic struc- 
ture calculations of solids and molecules. Its accuracy re- 
lies upon good approximations to the unknown exchange- 
correlation (xc ) energy as a functional of the electronic 
spin densities[2]- Over the past four decades, a vari- 
ety of increasingly sophisticated approximations has been 
developed [sj. The most commonly used in solid-state 
calculations today is the PBE version^] of the general- 
ized gradient approximation (GGA), employing both the 
density and its gradient at each point in space. Popular 
GGAs represent a well-tempered balance between com- 
putational efficiency, numerical accuracy, and reliability, 
but PBE also juggles the demands of quantum chemistry 
and solid-state physics[5i]. 

While PBE represented a high-point of non-empirical 
functional development 11 years ago, much has since 
been learned about its limitations. PBE reduces the 
chronic overbinding of the local spin density approxi- 
mation (LSDA)[H but, while LSD A often slightly un- 
derestimates equilibrium lattice constants by about 1%, 
PBE usually overestimates them by about the same 
amount. Other equilibrium properties, such as bulk mod- 
uli, phonon frequencies, magnetism, and ferro-electricity, 
are sensitive to the lattice constant, and so are also over- 
corrected by PBE0. Surface energies are too low in 
LSDA, but are made lower stiU by PBeQ. 

However, attempts to construct a better GGA face a 
Procrustean dilemma [sl: Those with an enhanced gra- 
dient dependence 0, [loj improve atomization and total 
energies, but worsen bond lengths, while more recent sug- 
gestions of a GGA for solids and surfaces 
|l5[ have a reduced gradient dependence and typically do 
improve lattice parameters and/or surface energies, but 
have been criticized for worsening energetics 13| . More 
advanced functionals have been constructed. For exam- 
ple, meta-GGA's, using also the orbital kinetic-energy 
densities, provide greater accuracy over a wider range of 
systems and properties [s!]. But current meta-GGA's do 



not improve lattice constants as dramatically as surface 
energies, and meta-GGA's are not yet available in all 
solid-state codes. 

In the present work, we explain the origin of this 
dilemma and show that, in a sense, no GGA can do 
both: Accurate atomic exchange energies require violat- 
ing the gradient expansion for slowly- varying densities, 
which is valid for solids and their surfaces. At the GGA 
level, one must choose. A pragmatic approach to lat- 
tice properties is therefore to use a modified functional 
especially for solids which, unlike previous suggestions, 
recovers the gradient expansion for exchange. Such a 
functional becomes exact for solids under intense com- 
pression. Analogous reasoning for correlation suggests 
fitting the jellium xc surface energy, as was first done by 
Armiento and Mattssonjl^. Our variation on PBE satis- 
fies these new conditions, by simply altering two param- 
eters in the original formula. We find our new functional 
significantly improves most lattice constants and surface 
energies while worsening atomization energies (as does 
AMOS of Ref. I4I). By restoring the gradient expansion, 
PBEsol also yields excellent jellium surface exchange en- 
ergies, and provides an improved starting point for more 
advanced functional construction. 

The GGA form for the exchange energy is simply 



d\eT\n{v))F^{s{v)) 



(1) 



where n(r) is the electronic density, e™'^(?T.) is the ex- 
change energy density of a uniform electron gas (propor- 
tional to n'*/^), s = \S/n\/{2kFn) (with kp = {'in'^nf/^) 
is the dimensionless density gradient, and Fx(s)_\s the 
so-called enhancement factor for the given GGA [J]. Eq. 
(1) is the spin-unpolarized form, from which the spin de- 
pendence can be deducedf^. Any GGA that recovers the 
uniform gas limit has 
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The gradient expansion 
varying electron gases has 



it! 



that is accurate for slowly- 
= 10/81 « .1234. 
To begin, Ref [17| showed that the exchange ener- 



gies of neutral atoms are very well approximated by 
their asymptotic expansion for large Z, i.e., ~ 
-.2208^5/3 „ iQQz + .... The first term arises from 
LSDA, but the second arises in a GGA from the con- 
tribution to Eq. (2) and requires jj, « 2^ge- Thus any 
GGA that is accurate for the exchange energies of neu- 
tral atoms must have /x « 2^qe- PBE does, although 
its value of /z = 0.219 was found from a different non- 
empirical argument. So does B88, as it was fitted to the 
X energies of noble gas atoms 18|. Even PW91 essentially 
does too [31, as only at irrelevantly small values of s does 
it revert to jicE- 

Thus, to attain accurate exchange energies of atoms 
(vital to dissociation energies in molecules and cohesive 
energies in solids), any GGA must strongly violate the 
gradient expansion for slowly- varying densities [l^. But 
most of thermochemistry occurs without free atoms, and 
is not much worse in LSDA than in PBE (e.g., [20|). 
Moreover, for the evaluation of exchange, the densities 
of real solids and their surfaces are often almost slowly- 
varying over space. Restoring the gradient expansion 
should improve their description (but worsen atomization 
energies). The GGA is a limited form, and cannot satisfy 
both conditions. Eq. (2) suggests a necessary condition 
for convergence of the second-order gradient expansion 
for exchange: s ^ 1. Since s < 1 for valence electrons 
in densely-packed solids (or s ^ 2 in core- valence regions 
of alkali atoms), and since the reduced Laplacian of the 
density is also < 1, the gradient expansion is important 
for exchange in solids. We choose ^ge for PBEsol. 

Now, for a GGA correlation functional that recovers 
the uniform gas limit, the gradient expansion is 
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d^rn(r){e™*(n(r))-H/3i"(r) + . . .} (3) 



where e™'^(n) is the correlation energy per particle of the 
uniform gas, /3 is a coefficient, and t = |Vn|/{2 fc-pF n} is 
the appropriate reduced density gradient for correlation 
(fixed by the Thomas-Fermi screening wavevector fcr f = 
\J Akp/T^, not kp.) For slowly- varying high densities [2l|, 
(icE = 0.066725. Unlike exchange, the second-order term 
in the gradient expansion for correlation cannot be small 
compared to the local term everywhere even for valence 
electrons in solids: fit^ can be large compared to |e™'^| 
(as PcEt^ = O.ls^/rs). The gradient expansion can be 
relevant to real systems (especially solids) for exchange, 
but much less so for correlation. 

The gradient expansion for correlation is less relevant 
to solids than is fxc{q) for the response of the uniform 
gas to a weak potential A cos(q • r). The exact fxc{q) is 
almost independent of q, up to 2kF [l^]- Thus LSDA, 
which produces a constant (the value at g = 0), yields 
an accurate approximation for q < 2kp. But any GGA 



with a non-zero xc contribution to second order in Vn 
produces a term quadratic in q. Since we are interested 
in weakly-varying valence electron densities in densely- 
packed solids, we wish to retain this excellent feature of 
LSDA. If 



(4) 



there is complete cancellation between beyond-LSDA 
X and c contributions, restoring LSDA response. 

In PBE, the gradient expansion for correlation is 
respected, i.e., (3 — Pge, and ^ « 2^ge satisfies 
Eq. (H]). This choice agrees well with the PW91 ex- 
change functional, and the numerical real-space cutoff 
construction [l^ of the GGA exchange energy, and yields 
highly accurate exchange energies of atoms. But we have 
already argued that /i « 2^qe is harmful for many con- 
densed matter applications. Once we choose ^ge for ex- 
change, we cannot recover simultaneously the GEA for 
correlation and the linear response of LSDA for a uni- 
form density. Exact satisfaction of Eq. (j4|) would yield 
(3 — 0.0375, but an increased value will satisfy another, 
more relevant constraint for solid-state applications. 

For correlation, large neutral jellium clusters are our 
paradigm, for which E^c — > e™'^y -I- e^™^A -|- .. as the 
radius grows, where e^™^ is the jellium surface xc energy, 
V the volume of the cluster and A its area. A GGA that 
recovers e^""^* will be correct in leading- and next-order 
for neutral jellium clusters as oo, in a similar way 

to modern exchange GGA's for neutral atoms. More- 
over, the surface energy is dominated by xc contributions 



and 



is really a bulk-like property, arising mainly 



(103% at Ts = 2) from a moderately-varying-density re- 
gion (with s < 1) inside the classical turning plane. 

We check that this condition is compatible with the 
restoration of the gradient expansion for exchange. Be- 
cause jellium clusters have a uniform bulk density and 
because most of the surface energy comes from within, 
the gradient expansion should be accurate. We find, at 
bulk density rg = 3, the errors of the surface exchange 
energy are: LSDA 27%, PBE -11% and PBEsol 2.7%. 

We next follow the lead of Ref. [H, and fit e|"/f to 
determine our correlation functional. But the jellium 
e|"''^ is not known exactly. Figure 1 shows surface en- 
ergy enhancements relative to LSDA. The likely "range 
of the possible" for e^'^'''' extends from TPSS meta-GGA 
[1, 0] or RPA-I- [2^ at the low end of what is possi- 
ble (in agreement with the most recent Quantum Monte 
Carlo calculations [3]) to the RPA-like Pitarke-Perdew 
(PP);25] value at the high end. The TPSS value prob- 
ably provides the best value that a GGA should try to 
achieve. We choose (3 = 0.046 and ^ ~ ^xge (within 
the PBE form) for PBEsol, to best fit the TPSS results. 
PBEsol should improve most surface energies over LSDA, 
whereas PBE worsens them. 

Thus we have violated Eq. @ in favor of good sur- 
face energies. But our value for [3 is considerably closer 
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FIG. 1: Ratio of calculated surface exchange-correlation to 
that of LSDA as a function of rs for various approximations. 

to that of the hnear response requirement (0.037) than 
that demanded by complete restoration of the gradient 
expansion (0.066). The hnear response of PBEsol is rea- 
sonably close to that of LSDA, and closer than that of 
the gradient expansion. 



density. But PBEsol also becomes exact for solids un- 
der intense compression, where real solids and their sur- 
faces become truly slowly-varying, and exchange domi- 
nates over correlation [TtJ. In Fig. 2, we plot the enhance- 
ment factors of PBE and PBEsol. For a spin-unpolarized 
(C = 0) density n = 3/(47rr|), we define Fj;c{rs,s) by 

E^c''^] = / dVer"(n(r))F.,(rs(r),.(r)) (5) 

The high-density (rs 0) limit is Fa,{s) of Eq. (1). The 
nonlocality or s-dependence of GGA exchange is dimin- 
ished from PBE to PBEsol, making the latter somewhat 
closer to LSDA. Over the whole range s < 1, the PBEsol 
Fx is close to H- ^ges^- The range < s < 3 is energeti- 
cally important for most properties of most real systems, 
while ^ s ^ 1 and 1 < rs < 10 are the ranges for 
valence-electron regions in many densely-packed solids. 



class 


LSDA PBE TPSS PBEsol 


n 

4 simple metals 
5 semiconductors 

5 ionic solids 
4 transition metals 
total 


lean error 
-9.0 2.9 5.3 -0.3 
-1.1 7.9 6.2 3.0 
-8.4 8.5 6.8 2.0 
-4.0 6.4 2.5 0.0 
-5.5 6.6 5.4 1.3 


mean 
4 simple metals 
5 semiconductors 

5 ionic solids 
4 transition metals 
total 


absolute error 
9.0 3.4 5.3 2.3 

1.3 7.9 6.2 3.0 

8.4 8.5 6.8 2.7 
4.0 6.4 2.7 1.9 
5.6 6.7 5.4 2.5 



TABLE 1: Errors in equilibrium lattice constants (in AxlO"'^) 
on our data set of 18 solids, relative to experiment with esti- 
mates of the zero-point anharmonic expansion removed [2^. 



FIG. 2: Enhancement factors of PBE and PBEsol, for spin- 
unpolarized systems, as a function of reduced density gradi- 
ent, for various values of rs. 

PBEsol is exact for the uniform gas, and highly accu- 
rate for large jellium clusters, both of moderate and high 



To test our functional, we employ a test set of 18 solids 
from Ref [2^. These come in four groups: Simple metals 
(Li,Na,K,Al), semiconductors (C,Si,SiC,Ge,GaAs), ionic 
solids (NaF,NaCl,LiCl,LiF,MgO), and transition metals 
(Cu,Rh,Pd,Ag). The set is not claimed to be representa- 
tive, but was chosen for the availability of basis functions 
and anharmonic corrections^. Our calculations use the 
Gaussian orbital periodic code of Ref [2^ , with basis sets 
of the same or higher quality. In Table I, we list both 
the mean errors and the mean absolute errors for lattice 
constants in LSDA, PBE, TPSS, and PBEsol. With the 
sole and marginal exception of SiC, LSDA makes lattice 
constants too short, as indicated by the negative mean 
errors. With the sole and marginal exception of Na, PBE 
makes lattice constants too long. Overall, the systematic 
PBE overestimate is close to the systematic LSDA under- 
estimate, as shown by the total mean absolute errors, and 
TPSS cures this very little. On the other hand, PBEsol 
greatly reduces this overestimate, by a factor of almost 4, 
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except for semiconductors, where LSDA is unsurpassed. 



error 


LSDA PBE TPSS PBEsol 


mean error 
mean abs. error 


3.35 0.54 0.18 1.56 
3.35 0.67 0.26 1.56 



TABLE IL Errors in atomization energies (eV) for the AE6 
set of molecules, using the 6-311+G(3df,2p) basis set. 



properties to be significantly improved over PBE, es- 
pecially for solids under compression. New pseudopo- 
tentials, compatible with PBEsol, must be generated to 
use PBEsol with pseudopotential codes. We thank NSF 
(CHE-0355405, CHE-0457030, and DMR-0501588) and 
OTKA for support. 



PBEsol is not expected to give good atomization ener- 
gies. In Table 2, we give the errors on the AE6 data set 
of molecules. These 6 molecules (SiH4, S2, SiO, C3H4 
(propyne), C2H2O2 (glyoxal), and C4H8 (cyclobutane) 
were chosen p 71] to be representative, i.e., to reproduce 
the errors of much larger data sets. As is clear, and ex- 
pected, PBEsol is much less accurate than PBE, only 
about halving the error of LSDA. This can be largely at- 
tributed to PBEsol's worsened total energies of atoms. 
For these atoms, the typical error in the total energy per 
electron is 0.2-0.3 eV for PBE, but 1.1-1.3 eV for PBEsol. 

We have demonstrated the relevance of the second- 
order gradient coefficient for the exchange energy of a 
slowly- varying density to the bulk and surface properties 
of sohds. The TPSS meta-GGA Q, which incorporates 
this coefhcient, gets good surface energies but its lattice 
constants are only marginally better than those of PBE 
on which it builds, whereas PBEsol is significantly better. 
This suggests that an improved meta-GGA might need to 
recover the gradient expansion for exchange over a wider 
range of n(r) than TPSS docs. 

Previous attempts to improve on PBE within the 
GGA form have retained the PBE gradient coefficients 
fi and 5 for small s, but altered the behavior at large 



s@, [^, 0, or have zeroed out /z 
fundamentally different from PBEsol. 
functional performs very similarly to 



11|,|12| 



and are thus 
The AM05[l3 
PBEsol for the 
solids studied here, but AMOS follows the proposal of 
Vitos et al.[ll| to fit the exchange energy density of an 
Airy gas (noninteracting electrons in a linear potential, 
LAG). This is not a model solid or surface density in any 
global sense. This approach does not recover the gradient 
expansion for exchange, because the energy density has 
no such expansion [28'|. For rg = 2 to 6, LAG has errors of 
11 to 56% for e|"'f , compared to 1.6 to 4.1% for PBEsol. 
By retaining the PBE formj^j, PBEsol also correctly re- 
tains the Lieb-Oxford bound and the high-density limit 
for correlation, while AMOS does not. Numerical com- 
parisons and details are available [2^. 

We have identified the simple exchange-correlation 
physics underlying many properties of many solids, and 
shown how it differs from that for atoms. We recommend 
PBEsol for the applications discussed here. Any exist- 
ing code that implements PBE can be instantly mod- 
ified to try PBEsol, by simply replacing the values of 
IJ, and /?. Modified PBE subroutines are available from 
pittp: / /dft . uci . edu . We expect geometries and related 
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